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ABSTRACT 

A simple and computationally efficient parameterization of the deleptonization, the entropy changes, 
and the neutrino stress is presented for numerical simulations of stellar core collapse. The param- 
eterization of the neutrino physics is based on a bounce profile of the electron fraction as it results 
from state-of-the-art collapse simulations with multi-group Boltzmann neutrino transport in spheri- 
cal symmetry. Two additional parameters include an average neutrino escape energy and a neutrino 
trapping density. The parameterized simulations reproduce the consequences of the delicate neutrino 
thermalization/diffusion process during the collapse phase and provide a by far more realistic alterna- 
tive to the adiabatic approximation, which has often been used in the investigation of the emission of 
gravitational waves, of multidimensional general relativistic effects, of the evolution of magnetic fields, 
or even of the nucleosynthesis in simulations of core collapse and bounce. For supernova codes that 
are specifically designed for the postbounce phase, the parameterization builds a convenient bridge 
between the point where the applicability of a stellar evolution code ends and the point where the 
postbounce evolution begins. 

Subject headings: supernovae: general — neutrinos — radiative transfer — hydrodynamics — methods: nu- 
merical 



1. INTRODUCTION 

Since the mid sixties, the gravitational collapse of 
stars has been studied base d on compute r simula- 
tions ijColgate fc White! 1196ft lArnetti Il96fift . Inves- 
tigated topics in the collapse phase inclu ded general 
relati vistic dynamics llMav fc White! Il967|) . magnetic 
fields l)Leblanc fc WilsorJl97flt). deleptonizat ion and neu- 
trino tra pping llSatol 119751 iMazurek 11976ft . progenitor 
rotati on iMuller fe HillebrandtJ 11981 ). dissociation en- 
ergy (IAxriejj[L9.g2j) 1 neutrino tra nsport and thermaliza- 
tion (Wilson 1985: lBr uennlll9§^ [l 1 nuclear electron car; 
ture ijBethe et alJl!979t iCooperstein fc Wambachlll984 " 



and the emission of gravitational waves ijMooreT 1981 
to name only a few of the many possible references. 
Later, the field has somewhat separated into two 
rather disjunct lines of research. On the one hand, 
the increasing computing power was focused on the 
details of the neutrino physics and neutrino trans- 
port in spherical symmetry. On the other hand, the 
computational resources were invested in multidimen- 
sional d ynamics for the investigation of rota ting pro- 
genitors JOtt et al.M2004 lArdelian et al]l2004l>. general 
relativity llDimmelmeier et alJ 12003 iSiebel et al.l l2003t 
iShibatafc Sekiguchill2005HDuez etal] 12005ft. and mag- 
netic fields llYamada fc Sawaill2004llKotake et "all 120041: 
ILiebendorfer et alJ 12005ft . Only few recent multidimen- 
sional collapse simulations made the effort to include 
neutrino physics. Thes e schemes are either very com- 
putationally expensive jBuras et all 12003 iMiiller et all 
120041: iWalder et al.ll2005|) or rely on simpli fications of the 
neutrino transport and its m icrophysics ijKotake et alJ 
20033 iFrver fc WarTeTIEiol . 

Among the neutrino physics that is most difficult 
to capture are accurate co mposition-dependent r ates of 
electron captures on nuclei l)Langanke et a l. 2003) in the 
early and medium phase of collapse and the competi- 



tion between neutrino diffusion and neu trino thermal- 
ization in the later pha se of the collapse l)Bruennlll985t 
IMvra fc BludmaiJll989ft . 

Electrons at high matter density are degenerate so 
that they are captured with large energies on bound 
or free protons. The produced high energy neutrinos 
first thermalize by electron scattering until their energy- 
dependent mean free path is large enough to make dif- 
fusion competitive. As the trapped neutrinos become 
degenerate themselves in the late stage of collapse, the 
ability to stream away at low energy actually becomes 
the bottleneck for further deleptonization. Most re- 
cent simulations of this thermalization-diffusion process 
are based on an individual treatment of different neu- 
trino energy bins with Boltzmann neutrino transport 
( Mez zaca ppa fc Bruennl Il993bt iBruenn fc Mezza cappal 
1 19971 iRampp fc JankaTteOOOl ILiebendorfer et all 120011: 
iThompson et alJ 120031: iHix et all 12003ft . The neutrino 
physics enters the equations of hydrodynamics in the 
form of three different source terms: as an electron frac- 
tion change rate, as an energy or entropy source, and as 
a source of acceleration by neutrino stress. 

This paper aims to bridge the two lines of research by 
an efficient and very simple prescription of how published 
and future results of accurate neutrino transport simula- 
tions in spherical symmetry could be incorporated in the 
hydrodynamics of core collapse for a more realistic study 
of the multidimensional dynamics, the role of magnetic 
fields, and the emission of gravitational waves than with 
adiabatic simulations. 

A parameterization of the three source terms is in- 
dividually described and evaluated in ^ 12141 Additional 
tests of the robustness of the parameterization with re- 
spect to model changes or differences in the dynamics 
are conducted in S|5] After the conclusion, more code- 
specific details of the implementation are defined in the 
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Fig . 1. — Electron fra ction profiles during core collapse in model 
G15 JEicbcndorfcr ct al. 2003). Each line shows the electron frac- 
tion as a function of density at a given time. The time slices have 
been chosen to represent each decade in the central density, p c , 
as indicated in the legend. The parameterization of the electron 
fraction, Y e , is based on the fact that the profile Y e (p) is only a 
weak function of time. 
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Fig. 2. — Comparison of the fit formula for the electron fraction 
profile at bounce with the original data of models N13 and G15. 
The agreement is much more accurate than the deviations inves- 
tigated in Fig. |3 The rise of Y e at the center of the G15 model 
is not reproduced by the fit. However, no improvement for the 
simulations would be obtained, because the minimum function in 
Eq. does not allow electron fraction increases anyway. 



appendix. 

2. DELEPTONIZATION 

Stellar core collapse proceeds by an imbalance between 
the self-gravitating forces of the inner core and its fluid 
pressure. The baryons contribute the most significant 
part to the gravitational mass of the stellar core while the 
degenerate electron gas provides the dominant contribu- 
tion to the pressure. The electron fraction, defined as the 
number of electrons per baryon in the gas, is therefore the 
most fundamental quantity for the stability of the inner 
core and the evolution of its size during the dynamical 
collapse. The electron fraction evolves by electron cap- 
tures on nuclei and the emission of the produ ced electron 
neutrinos. See ( Martin ez-Pinedo et al.ll2004j) for a recent 
review of the nuclear input physics before and during 
core collapse. The collapse continues until the matter at 
the center reaches nuclear density. Strong interactions 
reduce the compressibility at that point and the inner 
core bounces. The outgoing pressure wave turns into a 
shock wave as soon as it reaches the sonic point at the 
edge of the homologously collapsing inner core. The size 
of the inner core is important because it determines the 
location of this transition, the initial energy imparted to 
the shock, and the amount of matter outside of the shock 
that will be accreted and dissociated in the ongoing evo- 
lution. 

Figure ^ shows the electron fraction, Y^ 15 (p, t), as a 
function of density, p, at different times, t. The data has 
been taken from a general relativistic core collapse sim- 
ulation with Boltzmann neutrino transport and "stan- 
dard" input physics. The selected time slices correspond 
to the instances at which the central density reaches 10 10 
g/cm 3 , 10 11 g/cm 3 , . . ., 10 14 g/cm 3 , and finally p max at 
bounce. Figure^demonstrates that the function Y e (p, t) 
depends only weakly on time. Hence, it could be inter- 
esting to investigate how hydrodynamics simulations be- 
have when the computationally expensive calculation of 
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Parameters for the fitting-formula. 
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0.035 
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Y e (p, t) is replaced by linear interpolation in the logarith- 
mic density of a time-independent tabulated template 
of Y e (p). Because the electron fraction profile should 
be as accurate as possible at the time of bounce, when 
the final size of the inner core is determined, the choice 
Y e (p) = Y" e G15 (p, tt) at the time of bounce, tb, will be 
investigated. The data for the bounce electron fraction 
profile of model G15 are listed in machine-readable ta- 
bles in l)Liebendorfer et al.l20 05) . Alternatively, a fitting 
formula is provided here to increase the flexibility of the 
approach. The fitting of Y e (p) is based on a piecewise 
linear approximation with a piecewise cubic correction. 
The parameters are two points in the p-Y e space, (pi, Yi), 
and a scale, Y c , of the correction. Suggested values of the 
parameters are given in Table Q The fitting formula 
reads 

/ n r , . (, 2iogp-iogp 2 - log pi y 

x (p) = max — 1 . mm 1 , : : 

V logpa-logpi / 

Y e (x) = \(Y 2 + Y 1 ) + ^(Y 2 -Y 1 ) 

+ Y c [l-\x\+A\x\{\x\-l/2){\x\-l)]. (1) 

The comparison of the fit with the original data for mod- 
els N13 and G15 in Fig. 121 is very satisfactory. Tabled 
shows that the fits for the two models only differ in the 
density at the base of the silicon-oxygen layer, pi, and in 
the central electron fraction, Y 2 . The fit can thus easily 
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Fig. 3. — Electron fraction profiles as functions of enclosed mass 
at 10 ms and 2 ms before bounce, at bounce, and at 2 ms after 
bounce. The thick lines represent the accurate evolution of model 
G15 with Boltzmann neutrino transport. The thin lines represent 
the solution with the simple parameterization described in Eq. PI . 
Overall, the deleptonization is in agreement. The differences can be 
traced back to the earlier time slices where deviations of up to 5% in 
the deleptonization do occur. The time slice at 2 ms after bounce 
illustrates that the simple approximation for the deleptonization 
breaks down with the launch of the neutrino burst at an enclosed 
mass ~ 0.8 Mq. 



be modified by variation of these two physical quantities. 

The implementation of an electron fraction evolution 
along Y e (p) in a pure hydrodynamics scheme is achieved 
by the simple prescription 



6Ye = mm (0,Y e {p{t + St))- Y e (t)) 
St St 



(2) 



where the variation is taken at the same fluid element. 
The minimum function guarantees that the electron frac- 
tion monotonically d_ecreases even if transient instances 
would occur where Y e would be larger than the actual 
Y e . At the very start of a simulation, for example, Fig. 
n indicates that Y e < Y e almost everywhere (Y e in the 
progenitor data is represented by the profile belonging 
to the central density p c = 10 10 g/cm 3 ; Y e is represented 
by the profile belonging to p c — max). Thus, the delep- 
tonization described by Eq. (J2Jl sets in smoothly after 
a short time of adiabatic compression during which the 
original electron fraction profile moves to the right in the 
p-Y e space shown in Fig. \T\to join the bounce profile, Y e . 
The latter is then followed thereafter. 

Figure |3 compares the resulting evolution of the elec- 
tron fraction to the evolution obtained in the full model 
G15 with accurate neutrino transport. The small devi- 
ations between the two solutions are readily explained 
with the profiles shown in Fig. ^ At low density, early 
time slices (marked by triangles) have a lower electron 
fraction than the bounce profile (marked by a circle) . At 
densities ~ 10 11-12 g/cm 3 the realistic profile (marked 
by a diamond) assumes a larger electron fraction than 
the bounce profile. Time slices that reach to even larger 
densities (marked by a square and a cross) show lower 
electron fractions than the bounce profile. Because the 
bounce profile has been taken as template for the param- 
eterization, the same differences are reflected in Fig. [SJ 



most clearly visible in the profile at 2 ms before bounce. 
The parameterized deleptonization leads to somewhat 
higher Y e values in the outer layers, to lower Y e values 
in intermediate regions, and to higher Y e values near the 
center. The deviations are within 5%. 

The Y e value in model G15 rises again around nuclear 
density because of the increasing neutron degeneracy. 
Due to the minimum function, this change is not cap- 
tured by Eq. (J2J) so that the central Y e value at bounce 
in the parameterized evolution eventually falls below the 
corresponding value in the G15 model. A parameteriza- 
tion of the lepton fraction instead of the electron frac- 
tion would improve this, because the Y e would then as- 
sume its correct equilibrium value after neutrino trap- 
ping. Also the entropy changes would more consistently 
derive from lepton fraction changes than from electron 
fraction changes. The decision to work with the elec- 
tron fraction was motivated by the dominant role of Y e 
in the determination of the Chandrasekhar mass of the 
bouncing core, because Y e is a common input variable 
of realistic equations of state, and because a numerical 
determination of weak equilibrium would be required to 
find Y e otherwise. 

The Y e profile of model G15 at 2 ms after bounce 
shows the prominent electron fraction trough that arises 
when the accretion shock dissociates matter at neutrino- 
transparent densities around the enclosed mass ~ 0.9 
Mq . The huge number of neutrinos emitted in the neu- 
trino burst also leads to neutrino absorption ahead of 
the shock (manifest in the Y e peak at 0.95 M Q in the 
time slice at 2 ms after bounce). These Y" e -changes can 
of course not be expected to be represented by Eq. J5J 
based on the stationary bounce template. More sophisti- 
cated techniques are necessary to adequately implement 
the neutrino physics in the long-term postbounce phase. 

3. ENTROPY CHANGES 

The electron captures during collapse are not only 
changing the electron fraction, the matter entropy is af- 
fected as well. The baryons are in nuclear statistical 
equilibrium and the electrons are in thermal equilibrium. 
Thus, the increments of the entropy per baryon, Ss, are 
determined by the values of the chemical potentials p n , 
p p , and p e of neutrons, protons, and electrons, respec- 
tively. Additionally, there is an energy transf er between 
matter and neutrinos, Sq (e.g. l)Bruenn| [19 851 and refer- 
ences therein), 

TSs = -5Y e (fi e - p n + p p ) + Sq. (3) 

The temperature of the fluid is denoted by T. Depend- 
ing on the density of the material and the energy of the 
produced neutrinos, the neutrino can either (i) directly 
escape, (ii) thermalize and escape, or (hi) be trapped for 
longer than the dynamical time scale. 

In regime (i), Sq = SY e (E) where ( E) is the av erage 
energy of the freely escaping neutrinos l)Betheil990D . thus 

TSs = -SY e (p e -p n +p p - (E)) . (4) 

In this regime, electron capture on nuclei dominates over 
electron capture on prot ons. Due to th e average Q-value 
of the nuclei ~ 3 MeV l)Bruennlll985|) the neutrinos are 
produced with an average energy that is only marginally 
larger than p e —p n +p p . The corresponding small entropy 
decrease in this regime shall be neglected in the following 
parameterization. 
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Fig. 4. — Profiles of the entropy per baryon as functions of en- 
closed mass at 10 ms and 2 ms before bounce, at bounce, and at 2 
ms after bounce. The thick lines represent the accurate evolution 
of model G15 with Boltzmann neutrino transport. The thin lines 
represent the solution with the simple parameterization described 
in Eq. JFJ- The entropy evolution is in good agreement during 
collapse, but drops to about 7% below the reference values after 
bounce. The time slice at 2 ms after bounce shows again the promi- 
nent difference caused by the neutrino burst. The corresponding 
energy loss is not captured by the parameterization and leads to 
an overestimate of the entropy around an enclosed mass ~ 0.8 Mq. 



In regime (ii) the produced neutrinos start to be 
trapped by coherent scattering off heavy nuclei. The 
increasing matter density causes the electron chemical 
potential to rise so that the neutrinos are produced with 
higher initial energies. As the neutrino mean free path is 
proportional to E~ 2 , the fastest way of escape proceeds 
through thermalization at the high end of the neutrino 
energy spectrum with a tra nsition to diffusion and escape 
at its low end. Figure 7 in ijMartmez-Pinedo et alJl2004|) 
illustrates this thermalization process, which ends with a 
final escape energy of order Ef, sc ~ 10 MeV. For the fol- 
lowing parameterization of the entropy changes, Ef, sc is 
used as a constant parameter that defines where regime 
(ii) commences, namely when \i e — //„ + pt p — Ef, sc > 0. 
The entropy increase is then evaluated according to 



5s 
Hi 



5Y e /j, e - fi n + fj,p - E e v s 



(5) 



St T 

where 5Y e /6t is given by Eq. (J2J. 

At even higher densities, in regime (iii), the neutrinos 
are not able to escape before bounce. They are in equilib- 
rium with the fluid so that Sq in Eq. © is determined by 
the neutrino chemical potential, 5q = 5Y e p v . The neu- 
trinos form a normal gas component without transport 
abilities and the entropy is conserved in accordance with 
Eq. (3J and the equilibrium condition p e + p p = /i„ + \i v . 
As second parameter, a density threshold pt ra p ~ 2 x 10 12 
g/ cm 3 is introduced to define the beginning of regime (iii) 
beyond which no further entropy changes are taken into 
account. 

Figure 01 presents the entropy evolution based on Eq. 
(0 in comparison to the entropy profiles of model G15 
with comprehensive neutrino transport. The two pa- 
rameters used in Eq. (JSJ enable nice agreement in 



the collapse phase (profiles at 10 ms and 2 ms before 
bounce). The density at which the entropy starts to rise 
is mainly controlled by Eff c = 10 MeV. The level of the 
central entropy is mainly determined by the choice of 
Ptrap = 2 x 10 12 g/cm 3 . The values of the two parame- 
ters have not been particularly fine-tuned because they 
should rather represent a generic estimate than a multi- 
digit optimum for one particular model. The ~ 7% lower 
entropy in the parameterized evolution has no noticeable 
consequences for the dynamics of the cold matter in the 
inner core. The entropy profiles at 2 ms after bounce 
show the expected differences at late time due to the 
omission of the neutrino burst in the parameterization. 
The postshock region in the realistic G15 calculation is 
significantly cooler because of the energy loss inferred by 
the neutrino emission. The adjacent region ahead of the 
accretion shock is somewhat hotter in model G15 due to 
neutrino absorptions. 

4. NEUTRINO STRESS 

Neutrino stress is the third important quantity at the 
interface between neutrino transport and hydrodynam- 
ics. Although the neutrino pressure only contributes a 
fraction to the gas pressure, it influences the size of the 
inner core when the gas pressure and the gravitational 
forces cancel to a large extent. In regime (iii), i.e. the 
optically thick region where transport processes can be 
neglected, the neutrino stress is determined by the gra- 
dient of the neutrino pressure 



dv 
~~di 



P 



■ dpu 
dm ' 



(6) 



where dv/dt is the Lagrangean time derivative of the ve- 
locity. The term in the middle is the general expression; 
the term on the right hand side is its spherically sym- 
metric limit based on an enclosed mass m(r) at radius 
r. All general relativistic effects shall be neglected in the 
derivation of this simple parameterization. The neutrino 
pressure, p u , is readily evaluated based on the thermo- 
dynamic state of the fluid, 
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The neutrino chemical potential is given by p v = p, e — 
ji n + pp. The constants h, c, and k refer to the Planck 
constant, the speed of light, and the Boltzmann con- 
stant, respectively. F n (77) = / °° x n (e x ~ v + l)" 1 dx is 
the Fermi-Dirac function of order n l)Rhodeslll950() . 

For an implementation of neutrino stress in the opti- 
cally thin regime, an estimate of the neutrino number 
luminosity is needed. It can be derived from the delep- 
tonization in Eq. (J2J and the requirement of lepton con- 
servation. The simplest possible approximation in the 
construction of the neutrino flux from distributed sources 
is that the neutrinos leave isotropically and without time 
delay from the locations of deleptonization. Even if this 
assumption is obviously wrong in core collapse, it leads 
to a use ful first approximatio n of the non-local flux ge- 
ometry llGnedin fc Abell EoOl). Firstly, the deleptoniza- 
tion scheme in Eq. (JSJ already suppresses sources at 
high opacities. Secondly, lepton conservation requires 
that an isotropic source contributes with the square of 
its inverse distance so that closeby sources from this side 



5 



of the neutrinosphere influence the direction of the neu- 
trino flux more strongly than remote sources from the 
opposite side. With this assumption, the neutrino num- 
ber flux can be expressed by the gradient of a scalar field, 
ip, which fulfills the Poisson equation 



(8) 



In spherical symmetry, the neutrino number luminosity 
based on Eq. JHJ results in a straightforward integration 
of all enclosed sources, 



or 



i(r) 



^-N A dm. 
ot 



(9) 



Avogadro's number is denoted by Na- 

Let us stay in the spherically symmetric limit to de- 
rive an extension of the neutrino pressure to the optically 
thin regions. The neutrino stress can be expressed as a 
convolution of the neutrino number luminosity, the neu- 
trino en ergy, and the energy-dependen t reaction cross 
sections l|Mezzacappa fc Bruennlll993aj) . If the reaction 
cross sections are represented by the inverse mean free 
path, A -1 (E), one obtains based on the neutrino distri- 
bution function, / (E, p), 

Hii 9w r i 

-f(E,p)E 3 dEpdp. (10) 
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The neutrino momentum phase space element, E 2 dEdp, 
is described by the neutrino energy, E , and the propaga- 
tion angle cosine, /_t, with respect to the radial direction. 
Let's now assume that the spectrum of the number lu- 
minosity is well approximated by a Fermi-Dirac function 
with degeneracy parameter rj. The number luminosity 
can then be expressed as 



L = Anr c- 
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J f(E,p) E 2 dEpdp 



{kT„f F 2 (rj) H. 



(11) 



For the following rough estimate of the decline of the 
neutrino stress in the outer layers, it is further assumed 
that the mean free path scales as A^ 1 cx pE 2 , so that 
the energy integration in Eq. I|10|) can be evaluated and 
expressed by the number luminosity in Eq. jjlip. 

F 5 M L , . 

(12) 



The neutrino temperature, kT u , and the degeneracy pa- 
rameter, 77, cannot easily be derived as a function of ra- 
dius without solving a more detailed transport problem. 
Though, at high opacity, they limit to the matter tem- 
perature, kT, and the neutrino degeneracy, p v /kT, re- 
spectively. 

The neutrino stress for the hydrodynamics simulation 
is now generated by the following procedure: An estimate 
for the number luminosity is given by Eq. (PJ. Then, 
at all densities larger than a specified neutrino trapping 
density ptrap, the neutrino stress is evaluated according 
to Eqs. © and Q. With the neutrino stress at density 
Ptrap, a constant C is defined that attaches the scaling 
estimate found in Eq. (|12(l to the neutrino stress given 
byEq. ©, 



C = -Anr 



dm 



(kT) 



3 F 5 (p u /kT) L 



1 -1 



F 2 (pJkT) 4irr 2 



(13) 



In Eq. (|13fl . r, dp v /dm, T, p u , and L are evaluated 
at the transition density ptrap- Now, before Eq. Ijl2|l 
can be applied to extend the neutrino stress to the re- 
gion p < ptrap, approximations for kT v and 77 must be 
found for regime (i). There, it is assumed that the neu- 
trino spectrum is well represented based on a degeneracy 
r\ = 0. It is then consistent with 33 and the assumption 
of spherical symmetry if I express the neutrino temper- 
ature by the average neutrino energy parameter, E^ sc , 
and set kT v = F 2 (0) /F 3 (0) E e y sc . The simplest continu- 
ous transition of the neutrino stress from regime (iii) to 
regime (i) is realized by the adoption of the larger of the 
two limiting cases in each intermediate point. Hence, at 
p < Ptrap, the neutrino stress is approximated by 
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3 F 2 (0)F 5 (0) 



Fi(0) 

( 14 ) 

where r, T, p v = p e — p n + p p , and the estimate of 
L according to Eq. 10 represent the local values at the 
point where dv / dt is evaluated. Suggestions with respect 
to the implementation of this spherically symmetric neu- 
trino stress in a multidimensional hydrodynamics code 
and convenient approximations for the evaluation of the 
Fermi integrals are collected in the appendix. 

This procedure circumvents any explicit reference to 
cross sections and is therefore simple to apply. Of course, 
the cross sections do enter the constant C implicitly via 
the relation between the neutrino pressure gradient and 
the number luminosity at the transition density ptrap- As 
the value of the number luminosity is derived from Eq. 
(J2J), it reflects the deleptonization rate in the reference 
calculation with full transport which is based on infor- 
mation about the opacities used in the run that produced 
the electron fraction template at bounce time. 

Figure JSJ shows a summary of the evolution of veloc- 
ity profiles. It demonstrates that the collapse dynam- 
ics is accurately reproduced with the parameterized neu- 
trino physics. Some differences can still be made out: As 
discussed in ^ Eq. J5J implies a transient halt in the 
deleptonization until the density has increased enough 
that the initial Y e profile (labeled by p c = 10 10 in Fig. 
2]) catches up with the template (labeled by p c — max). 
This leads to a small delay in the infall of the outer layers 
with respect to the reference G15 model. This is visible 
in the velocity profiles at the enclosed mass of 1 - 1.2 
M Q . 

The time slice at 2 ms before bounce shows a soft out- 
going pressure wave in the solution with the parameter- 
ized neutrino physics that is not present in the reference 
model. The first suspicion was that it is related to the 
treatment of the neutrino stress because its appearance 
coincides to some extent with the moment when the den- 
sity ptrap is reached at the center, i.e. the moment when 
the applied neutrino stress jumps from zero to the value 
described by Eq. Q14[l. However, a more careful investi- 
gation showed that the dominant reason is the shallower 
electron fraction profile close to the center in the —2 ms 
time slice (see Fig. EJ. The outgoing wave is caused by 
a combination of the excess electron pressure at the cen- 
ter with the electron pressure deficit around an enclosed 
mass of 0.4 M Q . It introduces visible perturbations in 
the vicinity of the position of shock formation in the 
time slice at 0.2 ms before bounce. Nevertheless, the 
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Fig. 5. — Velocity profiles as functions of enclosed mass at 10 
ms, 2 ms, and 0.2 ms before bounce, at bounce, and at 2 ms after 
bounce. The thick lines represent the accurate evolution of model 
G15 with Boltzmann neutrino transport. The thin lines represent 
the s olut ion with the simple parameterization described in Eqs. JHJ 
and 1141 . The homologous infall is in nice agreement initially. The 
time slice at 2 ms before bounce shows a soft outgoing pressure 
wave that is mainly caused by electron fraction differences near 
the center. It causes perturbations around an enclosed mass of 0.6 
Mq in the —0.2 ms time slice. Otherwise, the strong pressure wave 
from bounce is in good agreement with the reference simulation at 
that time. The point of shock formation is very similar in both 
cases. The velocity profiles even continue to coincide some time 
beyond the emission of the neutrino burst as shown in the time 
slice at 2 ms after bounce. 



much stronger outgoing pressure wave from bounce and 
its steepening to the bounce-shock are nicely reproduced 
in the —0.2 ms and ms time slices. 

The significant differences induced by the neutrino 
burst in the electron fraction and entropy profiles at 2 
ms after bounce seem to have surprisingly little influence 
on the early dynamics. The velocity profiles as functions 
of enclosed mass still agree at 2 ms after bounce. The 
differences in the state of the postshock matter rather 
affect the shock radius than the shock mass. While the 
shock in the G15 model is positioned at a radius of 55 
km, it already has reached 63 km in the parameterized 
run due to the omission of the lepton and energy loss by 
the neutrino burst. 

5. MODEL DEPENDENCE 

After the demonstration that this simple parameter- 
ization works well for the one model G15, this section 
aims to explore to what extent the approximations are 
robust against variations in the dynamics or the initial 
conditions. A first test is the application of the same 
parameterization with the same parameters for E^ sc and 
Ptrav to a different stellar model. The G15 model dis- 
cusse d above is launched from the ( Wooslev fc Weaver! 
1995) model sl5s7b2. It h as a quite typical i ron core 
~1.3M , The run N13 in llLiebendorfer et al.ll2005D is 
based on the progenitor of ( No moto fc Hashimotdll98ft^ 
with an especially small iron core ~ 1.17 M Q . 

Figure compares different runs at the time of core- 
bounce. Figure 0a, shows the bounce profile Y e (p) — 



Y^ 13 (p, t),) of the N13 run that was used to specify the 
deleptonization according to Eq. (J2J) for the parameter- 
ized run A13. The bounce profiles of the N13 and A13 
runs are displayed in Fig. [^D-d. The quality of the ap- 
proximation is very similar to the one discussed above 
for the G15 model. The same choice of the parameters 
El sc and pt rap seems to fit also other stellar models. 

The dynamics in the N13 and G15 models is also dif- 
ferent. The former has been calculated with Newtonian 
hydrodynamics and 0(v/c) neutrino transport, the lat- 
ter with general relativistic dynamics and transport. As 
a second test I investigate how the parameterized solu- 
tion reacts if the bounce template for a Newtonian sim- 
ulation is taken from a general relativistic model of the 
same progenitor model. Thus, a bounce-template G13 
in Fig. 0i has been produced from a general relativistic 
calculation with neutrino transport. _A repetition of run 
A13 with the exchanged template, Y e (p) = Y^ 13 (p, tt) 
instead of Y e (p) — Yf 13 {p,h), leads to the results la- 
beled by X13G in Fig. HjjD-d. In spite of the significant 
dynamical differences in the runs that produced the tem- 
plates N13 (solid line) and G13 (dashed line) in Fig. HJi, 
the differences in velocity, electron fraction, and entropy 
profiles between runs A13 and X13G are barely visible. 
This means that the template for the parameterization 
of the microphysics can be extracted from a run whose 
dynamics differs from the run in which the parameteri- 
zation will be applied. 

A third test has been performed by the exchange of 
bounce templates between the progenitor models. First, 
a new template has been created with a Newtonian sim- 
ulation of the sl5s7b2 progenitor model labeled by N15 
in Fig. The parameterized Newtonian run X13P has 
then been launched from the 13 M Q progenitor model 
using this Y e (p) = Y™ 15 (p,tb) template. Figure 03- 
d shows a small deviation ~ 6% in the location of the 
shock formation. Although a given template from one 
progenitor still seems to be useful for many dynamical 
investigations that are launched from different progeni- 
tor models, it might be recommended to use a template 
that is derived from the same progenitor to achieve the 
best accuracy. A closer inspection of the differences re- 
veals the well-known fact that the final deleptonization 
in the inner core is essential for the location of the shock 
formation. Because the parameterization in Eq. J5J does 
not support electron fraction increases, it is the minimum 
Y e that counts. The minimum Y e in the templates of the 
N13 and G13 models in Fig. HJl coincide, while Y e in the 
N15 run reached slightly lower values. 

6. LIMITATIONS 

The described parameterization has been designed to 
provide a very efficient and reasonably accurate way to 
lead a hydrodynamics simulation from the onset of col- 
lapse of the supernova progenitor star to bounce. It has 
also been shown to be accurate for the initial rebound of 
the stellar core. However, the simplicity of the described 
approach entails several limitations. 

The signal of gravitational waves is not likely to be lim- 
ited to the short time interval around bounce that has 
traditionally been investigated l)Ott et alJ (|2004fl and ref- 
erences therein). In delayed, and especially in neutrino- 
driven explosions, the signal may decay after the first 
peak around bounce only to regain strength on a longer 
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Fig. 6. — Profiles at the time of bounce are shown to investigate the dependence of the parameterization on model variations, (a) Electron 
fraction templates for the N13, G13, and N15 model determined in simulations with full neutrino transport, (b-d) Velocity, electron fraction, 
and entropy profiles, respectively, for the parameterized simulations at bounce. The label N13 indicates the reference results with neutrino 
transport. The parameterized simulations A13, X13G, and X13P differ only by the Y" e -template that has been used. Template N13 has 
been used for run A13, G13 has been used for run X13G, and N15 has been used for run X13P. Swapping templates between Newtonian 
and general relativistic simulations introduced no significant differences, while a swapping of templates between progenitors lead to slightly 
different positions of the shock formation that correlate with the minimum Y e reached in the inner core. 



time scale of several hundreds of milliseconds, when fluid 
instabilities in the hot shocked matter between the pro- 
toneutron star and the shock front grow to asymme- 
tries on larger scale s ljThompsonl200ri[Scheck et al.l2004 
iMuller et al.ll200^) . The dynamics and the input physics 
for the description of the cold matter in the collapse 
phase are fundamentally different from the physics in the 
hydrostatic, but turbulent, hot mantle of the protoneu- 
tron star, which eventually determines the supernova ex- 
plosion after core-bounce. The suggested scheme has not 
been designed for the postbounce phase and does not de- 
scribe any of the weak interaction physics relevant after 
bounce. 

Firstly, the scheme handles only deleptonization and 
not neutrino transport. No attempt has been made to 
describe the very delicate neutrino energy deposition be- 



hind the stalled accretion shock. Neutrino heating at 50 
ms after bounce and beyond is thought to drive impor- 
tant fluid instabilities. Secondly, the scheme is only ap- 
plicable to situations where the deleptonization is caused 
by compression of cold matter. Even if the deleptoniza- 
tion in the condensing accreting material continues to 
be captured after bounce, major contributions from the 
emission of the neutrino burst at 2 ms after bounce and 
from neutrino diffusion in the dense core are missed by 
the simple parameterization. 

Therefore, if one extends an otherwise adiabatic code 
with the suggested parameterization for a simulation of 
the postbounce phase, one obtains the benefit that simu- 
lations pass through a significantly more realistic config- 
uration at bounce and that overly fast prompt explosions 
for small progenitors are avoided. But no improvement 



is obtained toward an adequate description of the impor- 
tant neutrino physics in the postbounce phase. 

Another question is the reliability of the scheme in the 
collapse phase when fast rotation rates are applied. The 
density distribution is then significantly less spherically 
symmetric. I think one can be somewhat optimistic be- 
cause the parameterization of the deleptonization during 
collapse in Eqs. J5J) and JSJ) relies more on the local 
density evolution than on the global geometry. This is 
supported by the low sensitivity of the parameterization 
to global dynamical changes investigated in ij5j A global 
asymmetry changes rather the spatial distribution of lo- 
cal conditions than the quality of the local conditions 
themselves. Relevant local conditions are: the time scale 
of compression, the density gradient, or the curvature of 
the isodensity surface at the point of investigation (as 
a measure of the deleptonizing volume per surface area). 
These quantities as function of the density may change by 
several percent with respect to a spherically symmetric 
scenario, but probabely not by orders of magnitude. If 
this speculation is true, the parameterization would still 
be approximately applicable. It is likely that accuracy 
is lost with increasing asphericity, but the improvement 
with respect to adiabatic simulations could still be sig- 
nificant. 

It would now be tempting to make also the treatment 
of the neutrino stress multidimensional. Equation (j^J 
already has a multidimensional form and Eq. I|14|) is 
readily generalized to 
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where ip is determined by Eq. (jHJ. Unfortunately, it is 
not clear, how C should be chosen in the multidimen- 
sional setting. Simple choices may imply discontinuities 
in the transition of the neutrino stress from Eq. JHJ to 
Eq. HI 5(1 that could introduce undesired artefacts in the 
velocity field. 

The spherically symmetric implementation of the neu- 
trino stress, as layed out in the appendix, is very simple 
and physically well-defined, but fails to represent any 
asymmetries in the shape or temperature of the neutri- 
nospheres and in the emitted neutrino luminosities. This 
might still be an acceptable limitation because the influ- 
ence of the neutrino stress (momentum deposition) on 
the collapse dynamics is not as crucial as the much dis- 
cussed neutrino heating (energy deposition) later in the 
postbounce phase. In the postbounce phase, which is 
not in the scope of this parameterization, asymmetries 
in the neutrino field can have a significant impact o n 
the dynamics of the shock revival ijKotake et al] l2003a'). 
On the other hand, stellar evolutio n models point to 
rather slowly rotating inner cores ffleger et all 120051 
IHirschi et alJ 1200^1 and corresponding collapse simula- 
tions with multidimensional neutrino transport have not 
show n evidence for sizab le asymmetries in the neutrino 
field llWalder et alJl2005l . 

7. CONCLUSION 

The dynamics of core collapse is dominated by elec- 
tron pressure. Dynamical simulations are only realistic 
if they can take account of electron captures on bound 
or free protons during infall. At increasing densities the 



electron captures get inhibited by neutrino phase space 
blocking so that the ability to thermalize and emit the 
produced neutrinos significantly contributes to the deter- 
mination of the final electron fraction in the inner core. 
The lower the electron fraction, the smaller is the mass of 
the core that bounces when nuclear densities are reached 
and the smaller is also the initial energy imparted to the 
outgoing shock. These well-known relationships ask for 
a careful inclusion of neutrino physics in simulations of 
stellar core collapse. The most accurate treatment of 
neutrin o transport has been developed in spherical sym- 
metry lIMezzaca ppa fc Messerl 119991 IRampp fc Jankal 
I2002T iThompson et all 120031: ILiebendorfer et all 120041: 
iSumivoshi et alJ l2005|) where it has been coupled to 
the most recent evaluation of electron capture rates 
on heavy nuclei lILanganke et all 120031: IHix et aT1l2003t 
iMarek et alJl2005j) . 

The previous sections present an embarassingly sim- 
ple and computationally efficient parameterization of the 
deleptonization in the collapse phase so that the conse- 
quences of past and future improvements of the neutrino 
physics should be straightforward to incorporate in mul- 
tidimensional simulations that focus on other aspects of 
core collapse. The scheme is based on electron fraction 
profiles as function of density at bounce, Y e (p), that have 
been produced by spherically symmetric simulations with 
full neutrino transport. A fitting formula is provided to 
conveniently represent these electron fraction profiles and 
to allow for adjustments to future developments. It is 
found that the time-evolution of the electron fraction in 
the simulations, Y e _^p,t), follows in reasonable approxi- 
mation the profile Y e (p). The corresponding parameter- 
ization brings the electron fraction at bounce, when it 
is most important for the shock dynamics, very close to 
the value in the original simulation. Unfortunately, it is 
not exactly equal because small differences in the time- 
dependence of the deleptonization will lead to density 
differences at bounce so that Y e (p) is not evaluated at 
exactly the same position as in the original simulation. 

Once the deleptonization is given, the entropy changes 
and neutrino stress are estimated based on two additional 
parameters, an average neutrino escape energy, Ef, sc — 
10 MeV, and a neutrino trapping density, pt rap = 2 x 10 12 
g/cm 3 . These parameters have not been fine-tuned, but 
may require adaption if the weak interaction physics is 
changed. 

In a comparison of the parameterized runs with the 
original ones it is found that the dynamics of core col- 
lapse is quite accurately reproduced. The parameter- 
ized neutrino physics presents a significant improvement 
with respect to adiabatic treatments. Its accuracy may 
even rival with neutrino transport approximations that 
neglect neutrino-electron scattering. 

However, some clearly visible inaccuracies do develop. 
The deleptonization of central zones at high densities 
proceeds slower than in the reference simulation. This 
leads to a weak outgoing pressure wave before bounce- 
density is reached. It propagates to the point of shock 
formation where it causes moderate perturbations. In- 
accuracies due to a deviation in the timing of the delep- 
tonization are also found in the hydrodynamic structure 
ahead of the shock. These deviations are a direct conse- 
quence of the differences between the function Y e (p) and 
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the time evolution of the electron fraction in the reference 
simulation. 

The evolution of the electron fraction and entropy is 
based on the local thermodynamic state of the fluid. This 
should allow an application also in models where the dy- 
namics moderately differs from the spherically symmet- 
ric case with which the electron fraction template has 
been produced. The less important neutrino stress has 
been approximated in a spherically averaged manner. An 
extension to handle asymmetries in the neutrino stress 
has been sketched, but its accuracy in models with es- 
sentially aspherical density distributions remains to be 
investigated. The sensitivity of the parameterized runs 
to changes in the model has been tested in two exper- 
iments that are accessible with comprehensive neutrino 
transport: Firstly, a template produced with general rel- 
ativistic dynamics has been used in a parameterized run 
with Newtonian gravity and, secondly, a template pro- 
duced with a 15 Mq progenitor star has been used in 
a parameterized run launched from a 13 Mq progenitor 
star with a very different structure of the iron core. 

It was found that the parameterization is not sensitive 
to the dynamical differences between general relativistic 
and Newtonian gravity. The choice of the progenitor for 
the template, however, has a small influence on the pa- 
rameterization. The extent of the deviation is directly 



determined by the difference between the central Y e val- 
ues in the templates. In the investigated case, the Y e 
difference translates to ~ 6% difference in the point of 
shock formation. I thus believe that one and the same 
template can be used with different progenitors for qual- 
itative studies of the dynamics. But a progenitor-specific 
template is recommended for more quantitative investi- 
gations. 

The parameterization has been designed for the col- 
lapse phase and not for the postbounce evolution. It 
does not contain any physics related to neutrino heat- 
ing. During collapse, it relies on the weak interaction 
physics used for the production of the bounce template. 
The nuclear and weak interaction physics during col- 
lapse is an interesting field of research. Progress in 
its adoption in spherically symmetric simulation s has 
been made l|Langanke et alJ 120031: iHix et alll2003f) and 
is likely to continue. Hopefully, this paper enables the 
efficient transfer of results from state-of-the-art neutrino 
transport simulations in spherical symmetry to simu- 
lations with more spatial degrees of freedom, where 
the implementation of comprehensive weak interaction 
physics together with the accurate solution of the energy- 
dependent multidimensional Boltzmann equation is yet 
computationally prohibitive. 



APPENDIX 
IMPLEMENTATION DETAILS 

The simulations with the parameterized neutrino physics are based on the same spherically symmetric general 
relativistic hydrodynamics code AGILE th at was used to solve the equations of hydrodynamics in the original pro- 
duction of the electron fraction templates ijLiebendorfer et al.l 2005). The interaction between neutrino physics and 
hydrodynamics proceeds through additional source terms in the conservation equations, i.e. YJr xt , e ext , and S ext in 
ijLiebendorfer et aT]l2002j) for electron fraction change rates, total energy change rates, an d momentum change rates , 
respectively. The simulations use the realistic Lattimer-Swesty equation of state version 2.7 l|Lattimer fc Swestvll99lJ) . 
Note that this version is likely to crash whenever one tries to enter the nuclear regime (eosflag=2) with a saved guess 
of the proton fraction obtained from the dissociated regime (eosflag=3). As convergence in the dissociated regime is 
robust, it is advisable to save guesses only when they are returned with eosflag=2. In the parameterized runs with 
AGILE, Eq. J5J is used to set Y^ xt . The heating rate e ext is found by numerical iteration of the equation of state 
until the entropy change rate specified by Eq. 10) is realized. The parameter E„ sa is set to 10 MeV. Note that the 
approximation Ss/5t = in the region where fi e — p n -I- p p — E^ sc < will still produce a rate e ext ^ when the electron 
fraction changes. At densities larger than pt ra p — 2 x 10 12 g/cm 3 , the spherically symmetric limit of Eq. © is used 
to calculate the neutrino stress S ext , while at lower densities, Eq. ifHI) is used. In the general relativistic simulations, 
the gravitational effect of neutrino energy and pressure has only been taken into account at densities p > ptrap- The 
neutrino pressure is evaluated by Eq. (JJJ) and the specific neutrino energy is set to e„ — 3p„/p. The gravitational 
effect of the neutrino luminosity was neglected. 

All time derivatives, d/dt, in Eqs. (|2I9I5I6I10I12I14|TB|) are Lagrangean, i.e. taken at the same mass element. Most 
hydrodynamics codes are discretized in space and not in mass; they rely on Eulerian time derivatives, d/dt. For the 
implementation of above equations in these schemes I suggest to use operator splitting. For example, the conservation 
equations of hydrodynamics should be straightforward to extend with a conservation equation for electron number, 

^ (pY e ) + V • (v P Y e ) = 0. (Al) 

In this first step, the electron fraction at location x is updated from Y e (x, i) to an intermediate value Y* (x, t + St) by 
the advection of electrons. The electron fraction update is completed in a second step by application of Eq. @ , 

SY e _ Y e (x, t + St)- Y e * (x, t + St) _ min (0, Y e (p (x, t + St)) - Y* (x, t + St)) 

~st ~ si " si • ( ' 

The entropy update can similarily be split into an entropy conserving hydrodynamics update (in whatever form it is 
realized in the hydrodynamics code) and a Lagrangean change of the specific entropy. For the latter, SY e /5t in Eq. 
© is substituted by the result of Eq. ljA2() . 

It is numerically stable to apply the neutrino stress in an operator split fashion as well because the neutrino pressure 
contributes only of order 10% to the total pressure. The simple spherical limit in Eq. JBJ is adequate if the deviations 
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from spherical symmetry arc small; the multidimensional form is advised otherwise. The extension of the neutrino 
stress to the optically thin regime is less straightforward because its derivation was based on some arguments that 
apply only in spherical symmetry. Eq. I|13[) . f°r example, extracts the opacities from the spherically symmetric run 
that produced the electron fraction templates. As long as the asphericity in the density distribution does not exceed 
half a density scale height, I recommend to use spherically averaged conditions of the multidimensional configuration 
to evaluate the neutrino stress. The integration of the de nsity over spheres results in an enclosed mass as function of 
the radius. The integration of the deleptonization in Eq. <A2(1 over spheres leads to a luminosity estimate according to 
Eq. © . Based on the spherically integrated density, energy density, and electron density the equation of state delivers 
the thermodynamic conditions used in Eq. (|14H to derive the neutrino stress for spheres with densities p < ptrap- 
The applicability of the multidimensional Eqs. JHJ and (|15fl in highly asymmetric situations was not numerically 
investigated because their assessment would require the corresponding reference simulations with multidimensional 
neutrino transport. In any case, it is important to abandon the application of the neutrino stress in optically thin 
regions after bounce before the shock reaches the density ptrap in order to prevent the growth of the constant C in 
Eq. (| 1 31) beyond limits. 

The Fermi inte grals in j)H are required f or degeneracy rj > 0. Convenient approximations to F2 (77) and F% (77) ha ve 
been taken from (EDStcin & Pethick llll98H) . while F 5 (77) has been derived along similar lines based on (Rhodes 1950): 

^ 2 (r,)~i(ry 3 +7r 2 r7)+|c(3)e- a " 

F 5 (T)) + 5vrV + W + ^f) ~ ^e"" (A3) 

with C (3) ~ 1.202 and a = 2vr 2 / (9C (3)) - 1 ~ 0.825. 
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